We will start off by reading our images and displaying them in full color.
library(EBImage)
library(CRImage)
setwd("/Volumes/SD/GoogleDrive/0.DA/Assignments/4")
image <- readImage('image1.jpg')
myImage <- readImage('cat.jpg')
display(image, method="raster")
display(myImage, method="raster")
Now we will convert the images to grayscale and display them.
grayImage <- EBImage::channel(image,"gray")
myGrayImage <- EBImage::channel(myImage,"gray")
#If there is masking of channel from other packages use the following one:
#grayImage <- EBImage::channel(image,"gray")
display(grayImage, method="raster")
display(myGrayImage, method="raster")
Compared to our colored image which had dimensions of 640 x 427 x 3, our new image has dimensions 640 x 427 x 1.
This is because for each pixel we had 3 values (RGB), but now we have only 1 value (grayscale).
We will convert our grayscale images to binary images using thresholds.
For this we will first create a function that does this.
#convert to Binary
threshold_values <- function(image, threshold){
image[image>threshold]=1
image[image<=threshold]=0
return (image)
}
bin02=threshold_values(grayImage,0.2)
bin04=threshold_values(grayImage,0.4)
bin06=threshold_values(grayImage,0.6)
bin08=threshold_values(grayImage,0.8)
bin09=threshold_values(grayImage,0.9)
par(mfrow = c(5,2))
plot(bin02)
title("Threshold = 0.2")
plot(image)
title(main = "Original")
plot(bin04)
title("Threshold = 0.4")
plot(image)
title(main = "Original")
plot(bin06)
title("Threshold = 0.6")
plot(image)
title(main = "Original")
plot(bin08)
title("Threshold = 0.7")
plot(image)
title(main = "Original")
plot(bin09)
title("Threshold = 0.9")
plot(image)
title(main = "Original")
bin02=threshold_values(myGrayImage,0.2)
bin04=threshold_values(myGrayImage,0.4)
bin06=threshold_values(myGrayImage,0.6)
bin08=threshold_values(myGrayImage,0.8)
bin09=threshold_values(myGrayImage,0.9)
par(mfrow = c(5,2))
plot(bin02)
title("Threshold = 0.2")
plot(myImage)
title(main = "Original")
plot(bin04)
title("Threshold = 0.4")
plot(myImage)
title(main = "Original")
plot(bin06)
title("Threshold = 0.6")
plot(myImage)
title(main = "Original")
plot(bin08)
title("Threshold = 0.7")
plot(myImage)
title(main = "Original")
plot(bin09)
title("Threshold = 0.9")
plot(myImage)
title(main = "Original")
Now we will apply some digital image processing functions.
We will first define a function that does contrast stretching:
contrast_stretching <- function(image,a,b,y_a,y_b,L, alpha, beta, gamma){
a = a / L
b = b / L
y_a = y_a / L
y_b = y_b / L
if(a > 1 || b > 1){
print("Error")
}
image[image < a] <- alpha * image[image < a]
image[image >= a & image < b ] <- beta * (image[image >= a & image < b ] - a) + y_a
image[image >= b] <- gamma * (image[image >= b] - b) + y_b
return(image)
}
Now we will apply this to our images:
contrast_stretched_img1 <- contrast_stretching(grayImage, a = 50, b = 150, y_a = 30, y_b = 200,L = 255, alpha=0.5, beta=5, gamma=20)
par(mfrow = c(2,2))
display(contrast_stretched_img1, method = "raster")
title("Contrasted")
display(grayImage, method = "raster")
title(main = "Original")
contrast_stretched_img1 <- contrast_stretching(myGrayImage, a = 50, b = 150, y_a = 30, y_b = 200,L = 255, alpha=0.5, beta=5, gamma=20)
par(mfrow = c(2,2))
display(contrast_stretched_img1, method = "raster")
title("Contrasted")
display(myGrayImage, method = "raster")
title(main = "Original")
We are first going to define the function that performs clipping:
clipping <- function(image,a,b,L,beta){
a = a / L
b = b / L
if(a > 1 || b > 1){
print("Error")
}
image[image < a] <- 0
image[image >= a & image < b ] <- beta * (image[image >= a & image < b ] - a)
image[image >= b] <- beta * (b - a )
return(image)
}
Now we can apply it to our images:
clipped_img1 <- clipping(grayImage, a = 0.4, b = 0.6, L = 1.0, beta=2)
par(mfrow = c(2,2))
display(clipped_img1, method = "raster")
title("Clipped")
display(grayImage, method = "raster")
title(main = "Original")
clipped_img1 <- clipping(myGrayImage, a = 0.4, b = 0.6, L = 1.0, beta=2)
par(mfrow = c(2,2))
display(clipped_img1, method = "raster")
title("Clipped")
display(myGrayImage, method = "raster")
title(main = "Original")
We first define the function that does compression:
compress <- function(image,c,L){
image <- c * log10(1 + image)
return(image)
}
Now we can apply it to our images:
compressed_img1 <- compress(grayImage, c = 50, L = 1.0)
par(mfrow = c(2,2))
display(compressed_img1, method = "raster")
title("Compressed")
display(grayImage, method = "raster")
title(main = "Original")
compressed_img1 <- compress(myGrayImage, c = 5, L = 1.0)
par(mfrow = c(2,2))
display(compressed_img1, method = "raster")
title("Compressed")
display(myGrayImage, method = "raster")
title(main = "Original")
We used a much lower value of c for our image because it was overall brighter. Very high values of c made the picture completely white otherwise.
We will now perform histogram equalization on our images.
First let’s plot the normal and equalized histogram for the original picture.
print('Image - histogram')
## [1] "Image - histogram"
equalized_grayImage = equalize(grayImage, range = c(0, 1), levels = 256)
par(mfrow = c(2,2))
hist(grayImage)
hist(equalized_grayImage)
plot(grayImage)
title("original")
plot(equalized_grayImage)
title("equalized")
Now we are going to see the same histograme for our image.
print('Our image - histogram')
## [1] "Our image - histogram"
equalized_grayImage = equalize(myGrayImage, range = c(0, 1), levels = 256)
par(mfrow = c(2,2))
hist(myGrayImage)
hist(equalized_grayImage)
plot(myGrayImage)
title("original")
plot(equalized_grayImage)
title("equalized")
We can notice that histogram equivalization lowered the brightness on our picture.
We will use a new picture for the this part.
image <- readImage('image2.jpg')
grayImage <- EBImage::channel(image,"gray")
display(x = image, method="raster")
We will visualize intensities across positions in a row of our matrix.
Then we will apply a filter to this and see how it looks.
First let’s start with the original picture:
row <- grayImage[200,]
par(mfrow = c(2,1))
plot(seq(length(row)),row, xlab='Position', ylab='Intensity', type="l")
title("Original")
filter9 <- matrix( c(1/9) ,nrow=1,ncol = 9)
row_filtered_9 <- filter2(x = row, filter = filter9)
plot(seq(length(row)),row_filtered_9, xlab='Position', ylab='Intensity', type="l")
title("Filtered")
par(mfrow = c(2,1))
row <- myGrayImage[200,]
plot(seq(length(row)),row, xlab='Position', ylab='Intensity', type="l")
title("Original")
filter9 <- matrix( c(1/9) ,nrow=1,ncol = 9)
row_filtered_9 <- filter2(x = row, filter = filter9)
plot(seq(length(row)),row_filtered_9, xlab='Position', ylab='Intensity', type="l")
title("Filtered")
We can clearly observe how filtering reduces the small variations and smoothes out the line.
We will apply a 2D blur filter that will average out the neighouring pixels.
par(mfrow = c(1,2))
filter2d <- matrix( c(1/81) ,nrow=9,ncol = 9)
grayImage_2d <- filter2(x = grayImage, filter = filter2d)
display(grayImage_2d, method = "raster")
title('Original image blurred')
filter2d <- matrix( c(1/81) ,nrow=9,ncol = 9)
grayImage_2d <- filter2(x = myGrayImage, filter = filter2d)
display(grayImage_2d, method = "raster")
title('Our image blurred')
Now we will apply a 1D Filter whose values sum to 1, and apply again the same filter but transposed to the convolved image from the first filter.
par(mfrow = c(1,2))
filter1d <- matrix( c(1/9), nrow=9, ncol = 1)
grayImage_1d10 <- filter2(x = grayImage, filter = filter1d)
display(grayImage_1d10, method = "raster")
title("1D filter on original image")
grayImage_1d10_transp <- filter2(x = grayImage_1d10, filter = t(filter1d))
display(grayImage_1d10_transp, method = "raster")
title("Extra 1D filter on previous (convolved) image")
We can notice that after the two filters the image is very similar to the image produced by the 2D filter.
But after doing some precise checks in R it doesn’t seem like they are the same.
par(mfrow = c(1,2))
filter1d <- matrix( c(1/9), nrow=1, ncol = 9)
grayImage_1d10 <- filter2(x = myGrayImage, filter = filter1d)
display(grayImage_1d10, method = "raster")
title("1D filter on original image")
grayImage_1d10_transp <- filter2(x = grayImage_1d10, filter = t(filter1d))
display(grayImage_1d10_transp, method = "raster")
title("Additional 1D filter on previous (convolved) image")
We will now apply a high pass filter to our data.
But after applying it, we will normalize the data back into a normal range (0 to 1).
par(mfrow = c(1,2))
hp_filter <- matrix(c(0,-1,0,-1,4,-1,0,-1,0),nrow=3,ncol=3)
hp_image <- filter2(grayImage, hp_filter)
min <- min(hp_image)
max <- max(hp_image)
hp_image <- (hp_image - min) / (max - min)
display(grayImage, method="raster")
title("Original image")
display(hp_image, method="raster")
title("High pass filtered image")
par(mfrow = c(1,2))
hp_filter <- matrix(c(0,-1,0,-1,4,-1,0,-1,0),nrow=3,ncol=3)
hp_image <- filter2(myGrayImage, hp_filter)
min <- min(hp_image)
max <- max(hp_image)
hp_image <- (hp_image - min) / (max - min)
display(myGrayImage, method="raster")
title("Original image")
display(hp_image, method="raster")
title("High pass filtered image")
The filter we applied is the presents the edges in the image or more precisely places where the (second derivative of) change of intensity is large.
We will use a new picture for the this part.
image <- readImage('image3.png')
grayImage <- EBImage::channel(image,"gray")
display(x = image, method="raster")
We will calculate the gradient magnitude after computing the local gradients.
To compute the local gradients we will use simple 1D arrays of (-1,0,1) and (1,0,-1) as shown on lecture slides.
filter_x <- matrix(c(-1,0,1),nrow=1,ncol=3)
filter_y <- matrix(c(1,0,-1),nrow=3,ncol=1)
gradient_x <- filter2(x = grayImage, filter = filter_x)
gradient_y <- filter2(x = grayImage, filter = filter_y)
gradient_mag <- sqrt(gradient_x^2+gradient_y^2)
#display(x = gradient_x, method="raster")
#display(x = gradient_y, method="raster")
print("Gradient magnitudes of Lena")
## [1] "Gradient magnitudes of Lena"
display(x = gradient_mag ,method="raster")
title("Lena")
Now let’s explore the gradient magnitude at various thresholds:
First we start with Lena:
grad_mag_1 <- threshold_values(gradient_mag,0.05)
grad_mag_2 <- threshold_values(gradient_mag,0.25)
grad_mag_3 <- threshold_values(gradient_mag,0.5)
grad_mag_4 <- threshold_values(gradient_mag,0.75)
grad_mag_5 <- threshold_values(gradient_mag,1.0)
par(mfrow = c(5,2))
plot(grad_mag_1)
title("Threshold = 0.05")
plot(grayImage)
title(main = "Original")
plot(grad_mag_2)
title("Threshold = 0.25")
plot(grayImage)
title(main = "Original")
plot(grad_mag_3)
title("Threshold = 0.5")
plot(grayImage)
title(main = "Original")
plot(grad_mag_4)
title("Threshold = 0.75")
plot(grayImage)
title(main = "Original")
plot(grad_mag_5)
title("Threshold = 1.0")
plot(grayImage)
title(main = "Original")
We can notice that only stronger (intensive) edges show up when we apply higher thresholds.
Now we will do the same for our image:
filter_x <- matrix(c(-1,0,1),nrow=1,ncol=3)
filter_y <- matrix(c(1,0,-1),nrow=3,ncol=1)
gradient_x <- filter2(x = myGrayImage, filter = filter_x)
gradient_y <- filter2(x = myGrayImage, filter = filter_y)
gradient_mag <- sqrt(gradient_x^2+gradient_y^2)
display(x = gradient_mag ,method="raster")
title("My image")
grad_mag_1 <- threshold_values(gradient_mag,0.05)
grad_mag_2 <- threshold_values(gradient_mag,0.25)
grad_mag_3 <- threshold_values(gradient_mag,0.5)
grad_mag_4 <- threshold_values(gradient_mag,0.75)
grad_mag_5 <- threshold_values(gradient_mag,1.0)
par(mfrow = c(5,2))
plot(grad_mag_1)
title("Threshold = 0.05")
plot(myGrayImage)
title(main = "Original")
plot(grad_mag_2)
title("Threshold = 0.25")
plot(myGrayImage)
title(main = "Original")
plot(grad_mag_3)
title("Threshold = 0.5")
plot(myGrayImage)
title(main = "Original")
plot(grad_mag_4)
title("Threshold = 0.75")
plot(myGrayImage)
title(main = "Original")
plot(grad_mag_5)
title("Threshold = 1.0")
plot(myGrayImage)
title(main = "Original")
The value of 0.25 as a threshold seems rather good for both images.
Higher thresholds are not very informative because gradient maginuted rarely have values above those thresholds.
Now we will implement a function that allows the user to choose between Prewitt and Sobel operators, and allows the user to specify thresholds, as well as orientation (horizontal, vertical or both).
edge_detect <- function(image, type, orientation, threshold = NULL){
matrix_prewitt_x = matrix(c(-1, -1, -1, 0, 0, 0, 1, 1, 1), nrow = 3, ncol = 3)
matrix_prewitt_y = t(matrix_prewitt_x)
matrix_sobel_x = matrix(c(-1, -2, -1, 0, 0, 0, 1, 2, 1), nrow = 3, ncol = 3)
matrix_sobel_y = t(matrix_sobel_x)
filter_x = 0
filter_y = 0
if(type == "prewitt"){
filter_x = matrix_prewitt_x
filter_y = matrix_prewitt_y
}
else if (type == "sobel"){
filter_x = matrix_prewitt_x
filter_y = matrix_prewitt_y
}
else{
print("Error")
}
if(orientation == "both"){
gradient_x <- filter2(x = image, filter = filter_x)
gradient_y <- filter2(x = image, filter = filter_y)
gradient_mag <- sqrt(gradient_x^2+gradient_y^2)
}
else if (orientation == "vertical"){
gradient_y <- filter2(x = image, filter = filter_y)
gradient_mag <- sqrt(gradient_y^2)
}
else if (orientation == "horizontal"){
gradient_x <- filter2(x = image, filter = filter_x)
gradient_mag <- sqrt(gradient_x^2)
}
else{
print("Error")
}
if(! is.null(threshold)){
gradient_mag <- threshold_values(gradient_mag,threshold)
}
return (gradient_mag)
}
Let’s test the funtion now.
Sobel operator with both orientations and no threshold.
par(mfrow = c(1,2))
sobel_image <- edge_detect(grayImage, type = "sobel", orientation = "both", threshold = NULL)
display(sobel_image,method="raster")
sobel_image <- edge_detect(myGrayImage, type = "sobel", orientation = "both", threshold = NULL)
display(sobel_image,method="raster")
Sobel operator with horizontal orientation and no threshold.
par(mfrow = c(1,2))
sobel_image <- edge_detect(grayImage, type = "sobel", orientation = "horizontal", threshold = NULL)
display(sobel_image,method="raster")
sobel_image <- edge_detect(myGrayImage, type = "sobel", orientation = "horizontal", threshold = NULL)
display(sobel_image,method="raster")
Sobel operator with vertical orientation and no threshold.
par(mfrow = c(1,2))
sobel_image <- edge_detect(grayImage, type = "sobel", orientation = "vertical", threshold = NULL)
display(sobel_image,method="raster")
sobel_image <- edge_detect(myGrayImage, type = "sobel", orientation = "vertical", threshold = NULL)
display(sobel_image,method="raster")
Sobel operator with both orientations and threshold of 0.25.
par(mfrow = c(1,2))
sobel_image <- edge_detect(grayImage, type = "sobel", orientation = "both", threshold = 0.25)
display(sobel_image,method="raster")
sobel_image <- edge_detect(myGrayImage, type = "sobel", orientation = "both", threshold = 0.25)
display(sobel_image,method="raster")
Sobel operator with both orientations and threshold of 0.5.
par(mfrow = c(1,2))
sobel_image <- edge_detect(grayImage, type = "sobel", orientation = "both", threshold = 0.5)
display(sobel_image,method="raster")
sobel_image <- edge_detect(myGrayImage, type = "sobel", orientation = "both", threshold = 0.5)
display(sobel_image,method="raster")
Sobel operator with both orientations and threshold of 0.75.
par(mfrow = c(1,2))
sobel_image <- edge_detect(grayImage, type = "sobel", orientation = "both", threshold = 0.75)
display(sobel_image,method="raster")
sobel_image <- edge_detect(myGrayImage, type = "sobel", orientation = "both", threshold = 0.75)
display(sobel_image,method="raster")
Prewitt operator with both orientations and no threshold.
par(mfrow = c(1,2))
prewitt_image <- edge_detect(grayImage, type = "prewitt", orientation = "both", threshold = NULL)
display(prewitt_image,method="raster")
prewitt_image <- edge_detect(myGrayImage, type = "prewitt", orientation = "both", threshold = NULL)
display(prewitt_image,method="raster")
Prewitt operator with vertical orientations and no threshold.
par(mfrow = c(1,2))
prewitt_image <- edge_detect(grayImage, type = "prewitt", orientation = "vertical", threshold = NULL)
display(prewitt_image,method="raster")
prewitt_image <- edge_detect(myGrayImage, type = "prewitt", orientation = "vertical", threshold = NULL)
display(prewitt_image,method="raster")
Prewitt operator with horizontal orientations and no threshold.
par(mfrow = c(1,2))
prewitt_image <- edge_detect(grayImage, type = "prewitt", orientation = "horizontal", threshold = NULL)
display(prewitt_image,method="raster")
prewitt_image <- edge_detect(myGrayImage, type = "prewitt", orientation = "horizontal", threshold = NULL)
display(prewitt_image,method="raster")
Prewitt operator with both orientations and threshold of 0.25.
par(mfrow = c(1,2))
prewitt_image <- edge_detect(grayImage, type = "prewitt", orientation = "both", threshold = 0.25)
display(prewitt_image,method="raster")
prewitt_image <- edge_detect(myGrayImage, type = "prewitt", orientation = "both", threshold = 0.25)
display(prewitt_image,method="raster")
Prewitt operator with both orientations and threshold of 0.5.
par(mfrow = c(1,2))
prewitt_image <- edge_detect(grayImage, type = "prewitt", orientation = "both", threshold = 0.5)
display(prewitt_image,method="raster")
prewitt_image <- edge_detect(myGrayImage, type = "prewitt", orientation = "both", threshold = 0.5)
display(prewitt_image,method="raster")
Prewitt operator with both orientations and threshold of 0.75.
par(mfrow = c(1,2))
prewitt_image <- edge_detect(grayImage, type = "prewitt", orientation = "both", threshold = 0.75)
display(prewitt_image,method="raster")
prewitt_image <- edge_detect(myGrayImage, type = "prewitt", orientation = "both", threshold = 0.75)
display(prewitt_image,method="raster")
We can notice both the Prewitt and Sobel operators do a very good job.
I would even say they performed better than the gradient magnitude approach because with higher threshold some values still existed.
We will use the Laplacian operator to calculae the second derivatives.
I will try a mask in the form: [0,1,0; 1,0,1; 0,1,0] that could be used to detect zero crossings.
zero_crossing = function(image){
laplace_operator <- matrix(c(0, 1, 0, 1, -4, 1, 0, 1, 0), nrow = 3, ncol = 3)
filter_operator <- matrix(c(0, 1, 0, 1, 0, 1, 0, 1, 0), nrow = 3, ncol = 3)
zero_crossings_raw <- filter2(x = image, filter = laplace_operator)
zero_crossing_raw2 <- filter2(x = zero_crossings_raw, filter = filter_operator)
return(zero_crossing_raw2)
}
We can apply this function now to our images.
par(mfrow = c(1,2))
display(grayImage,method="raster")
display(zero_crossing(grayImage),method="raster")
par(mfrow = c(1,2))
display(myGrayImage,method="raster")
display(zero_crossing(myGrayImage),method="raster")
The results aren’t that good.
But we could improve the results using a Laplacian of Gaussian where Gaussian filtering is used to smooth out the transition and the Laplacian operator is used for edge detection.
We will use the same mask.
zero_crossing = function(image){
laplace_operator <- matrix(c(
0,0,-1,0,0,0,-1,-2,-1,0,-1,-2,16,-2,-1,0,-1,-2,-1,0,0,0,-1,0,0
), nrow = 5, ncol = 5)
filter_operator <- matrix(c(0, 1, 0, 1, 0, 1, 0, 1, 0), nrow = 3, ncol = 3)
zero_crossings_raw <- filter2(x = image, filter = laplace_operator)
zero_crossing_raw2 <- filter2(x = zero_crossings_raw, filter = filter_operator)
return(zero_crossing_raw2)
}
We can apply this function now to our images.
par(mfrow = c(1,2))
display(grayImage,method="raster")
display(zero_crossing(grayImage),method="raster")
par(mfrow = c(1,2))
display(myGrayImage,method="raster")
display(zero_crossing(myGrayImage),method="raster")
I think this approach with Laplacian of Gaussian gave pretty good results.